Coarse-grained interaction potentials for polyaromatic hydrocarbons 
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Using Kohn-Sham density functional theory (KS-DFT), we have studied the interaction between 
various polyaromatic hydrocarbon molecules. The systems range from mono-cyclic benzene up 
to hexabenzocoronene (hbc). For several conventional exchange-correlation functional^ potential 
energy curves of interaction of the 7T-7T stacking hbc dimer are reported. It is found that all pure local 
density or generalized gradient approximated functionals yield qualitatively incorrect predictions 
regarding structure and interaction. Inclusion of a non-local, atom-centered correction to the KS- 
Hamiltonian enables quantitative predictions. The computed potential energy surfaces of interaction 
yield parameters for a coarse-grained potential, which can be employed to study discotic liquid- 
crystalline mesophases of derived polyaromatic macromolecules. 



I. INTRODUCTION 

Discotic thermotropic liquid crystals can be formed by 
flat molecules with a central aromatic core and several 
aliphatic chains attached at the edges 0,0- The size and 
the shape of the cores can be varied, as well as the length 
and the structure of the side chains, which allows the 
control of functional properties of these mesophases 0, 
0, 0. Liquid-crystalline properties such as fluidity are 
of help to process these compounds and even to develop 
self-healing materials. 

The self-organization of the aromatic cores into 7r-7r 
electron bonded stacks surrounded by saturated hydro- 
carbons allows one-dimensional charge transport along 
the columns 0, 0] • Unfortunately, the spatial arrange- 
ment of stacks is not perfect, i.e. the columns can be 
misaligned, tilted, or form various types of topological 
defects. In addition, the local alignment of molecules 
in columns can vary for different compounds. This con- 
siderably affects the intermolecular overlap of the it or- 
bitals and thereby the efficiency of the charge transport 
in a single column. As a consequence, the details of the 
morphology of the conducting film are crucial 0, 0, ^| . 
An accurate in silico prediction of mesophase properties 
prior to the actual synthesis of the compounds could ac- 
count for a considerable gain in efficiency on the route 
towards the design of macromolecular photo devices • 

An accurate understanding of the constituting 
molecules and their intermolecular interactions is manda- 
tory for controling the local alignment of the disks or the 
global arrangement of the columns in the mesophase. De- 
pending on the chosen length- and time-scales different 
methods can be considered. In principle, at the quan- 



tum chemistry level, one can study electronic, inter and 
intra molecular adsorption and adhesion processes |12| . 
and even compound design without empiricism. Still 
at an atomistic resolution - but using empirical molecu- 
lar force-fields - nanometer and nanosecond simulations 
can yield local properties, such as order parameters, or 
molecular arrangements [La. fl5L Ha . fl7| . In an even more 
extended time and length scale limit (/im and /us), coarse- 
grained simulations allow to describe the morphology of 
bulk material, global arrangement of macromolecular ob- 
jects such as columns or generic phase diagrams, and de- 
fects mHumminiii. 

To our knowledge, discotic materials have been stud- 
ied only very little and if so with idealized model poten- 
tials. The reason being that an ah initio treatment of 
the dispersion forces is computationally very demanding. 
Molecular dynamics (MD) is able to treat larger systems, 
however the details of the electronic structure, which are 
crucial for an understanding of electron transport, are by 
construction not included in MD. Moreover, it requires 
empirical parameters for the atom-atom interactions, and 
defects and the mesophase morphology can only be stud- 
ied at even more coarse levels. Consequently, multiscale 
methodologies [U |H EHH HI M OHl seem to repre- 
sent the most adequate and tractable description of these 
systems. 

The aim of this work is to make first steps towards 
multiscale modeling of discotic mesophases of polyaro- 
matic hydrocarbons. Namely, coarse-grained potentials 
for the interaction between representative polyaromatic 
molecules in a face-to-face geometry are obtained from 
first principles. They can be applied to the study of 
macroscopic properties of these materials or their chem- 
ically derived structures. 



II. COMPUTATIONAL DETAILS 
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These forces result from the correlated fluctuation of 
non-overlapping electron densities of molecular frag- 
ments Their prediction from first principles has 
remained a long-standing challenge because of the very 
high accuracy required to describe electron correlation 
effects. Explicitly correlated wavefunction methods such 
as coupled-cluster, configuration interaction, or quan- 
tum Monte Carlo allow for an accurate treatment of 
these forces but are computationally prohibitively ex- 
pensive for all but the smallest polyaromati c hy drocar- 
bons, such as for instance the benzene dimer J33|. Kohn- 
Sham density functional theory (KS-DFT), on the other 
hand, would be an exact electronic structure method 
if the true exchange-correlation (xc) term in the KS- 
potential was known. Unfortunately, this is not the case 
and for all practical purposes approximations have to 
be made. While some of the pure xc functionals for- 
tuitously but inconsistently predict binding for London- 
dispersion complexes, it is not yet generally possible to 
describe correctly vdW interactions within DFT using 
the local density approximation (LDA), the generalized 
gradient approximation (GGA) or even the - on aver- 
age more accurate - hybrid exchange-correlation func- 
tionals nnummiHiiii. 

Considering the ubiquitous nature of these intermolec- 
ular forces and their importance for self-assembly, much 
effort is being devoted to design superior xc-potcntials 
which can account correctly for all intermolecular inter- 
actions. The use of nonlocal correlations by electron den- 
sity partitioning |4Ct |41| can efficiently remedy this de- 
ficiency but implies an a priori assignment of molecular 
fragments. A 'van der Waals' functional as proposed in 

Ref. nils nam, and based on response theory be- 
comes rapidly intractable, such as the schemes described 
in Ref. [471 l4q . As a consequence, empirical a posteriori 
pairwise atom-atom based correction terms J37l to 
the energy are in wide spread use for practical applica- 
tions. The required parameters and damping functions 
for the correct repulsive behavior can be obtained from 
experiments or from various theoretical approaches in- 
cluding time dependent DFT [5(J, I^H 02, 03 • However, 
these r _6 -dependent corrections to the energy and ionic 
forces need artificial damping functions to allow for the 
correct repulsive behavior, and more importantly leave 
the electronic structure uncorrected. 

In this study, London-dispersion forces are computed 
from an improved electronic structure calculation which 
exploits a recently presented semi-empirical dispersion 
calibrated atom-centered (DCACP) correction, vf lsp , to 
a given Hamiltonian j54[. Specifically, it can be seen as 
a nonlocal extension of a given, local for LDA or GGA, 
xc-potential, 

v? c tended (r) = Mr) + £«f sp (r,r',{A}), (1) 

i 

where index i enumerates all atoms. As generally sug- 
gested in Ref. j55| . the atom-type-dependent parameters 
{A} require preliminary calibration for an improved elec- 




FIG. 1: Sketches of the polycyclic aromatic hydrocarbons 
benzene (A), pyrene (B), triphenylene (C), perylene (D), 
coronene (E), and hexabenzocoronene (hbc) (F). All hydro- 
gens are omitted for the sake of clarity. Bonds are represented 
by edges, carbon atoms by vertices. 



tronic structure fulfilling additional requirements, such as 
exerting a London-dispersion force on the ions. Specifi- 
cally, the BLYP-DCACP for carbon as it has been intro- 
duced, calibrated to the benzene dimer, and assessed in 
Refs. [54ll5rll57j is employed. Generalization of this cor- 
rection to other xc- functionals than BLYP has already 
been carried out jS^ • 

All DFT calculations have been carried out using 
the plane-wave basis set electronic structure pro gram 
CPMD 3.92 [13, the xc-functionals BLYP [6(J M H2, 
BP i3L PBE IH, LDA (using the Perdew and 
Zunger fit |65j to the data of Ceperley and Alder [6tijp. 
Goedecker pseudopotentials from Refs. |67l l68l Ir39| . and a 
plane-wave cutoff of 100 Ry. The isolated system module 
in CPMD has been employed together with the Poisson- 
solver of Tuckerman and Martyna F70| . The box-size is 
sufficiently converged at 19 x 19x20 A 3 for the largest sys- 
tem (hbc) and has been kept fixed for all molecules and 
all distances. Carbon-type DCACPs ^54| have been used 
only as a correction to the BLYP functional, no correc- 
tion has been employed for hydro gen atoms. For the cal- 
ibration of the DCACP's in Ref. [Hi], the M0ller-Plesset 
second order perturbation theory energy of interaction 
of the benzene dimer was used as a reference. Since the 
plane wave basis is independent of atomic positions, no 
basis set superposition errors occurr. Relative energies 
have been computed for identical box-sizes and cutoffs. 
For all geometry optimizations the residual tolerance for 
ionic forces has been set to 0.0005 a.u. For all calcula- 
tions of energies of interaction, the monomer geometries 
have been optimized with the given xc-functional. There- 
after, for the calculation of the interaction curves, the in- 
tramolecular geometries of the top-on-top moieties have 
been hold fixed and only the intermolecular distance has 
been varied. 

Several aromatic disks with symmetric cores have been 
selected. Sketches of their chemical composition are 
shown in Fig. 2] Namely, complexes of benzene (A), 
pyrene (B), triphenylene (C), perylene (D), coronene (E), 
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FIG. 2: Total potential energy of interaction for hexabenzo- 
coronene (see structure F, Fig.0 as a function of intermolec- 
ular distance for the LDA, BLYP, PBE, BP, and BLYP + 
DCACP xc-functionals. 



and hexabenzocoronene (hbc) (F) have been studied, 
which all fulfill Huckel's (4JV + 2) rule. 

III. RESULTS AND DISCUSSION 

A. Different density functionals 

First, several xc functionals have been assessed by eval- 
uating the potential energies of interaction, 

^int ^dimcr 2^ monomcr 

for hexabenzocoronene using LDA, BLYP, PBE, and BP. 
The results are shown in Fig. [3 As expected, LDA 
exhibits some fortuitous binding for graphite-like struc- 
tures, however when using the more sophisticated GGA- 
functionals, a completely repulsive behavior is obtained. 
It is found that the PBE gives the least repulsive in- 
teraction, followed by BP, and BLYP, suggesting that 
Becke's exchange potential is too repulsive for the in- 
vestigated systems. Upon inclusion of the DCACP- 
extension into the BLYP-DFT Hamiltonian, the interac- 
tion energy curve is in a reasonable agreement with what 
can be expected from experimental results, available for 
coronene [71] . 

B. Interaction energy profiles 

The interaction energy profiles E mt (r) have been cal- 
culated for all systems in the face-to-face geometry using 
the BLYP-DCACP KS-Hamiltonian. The results are pre- 
sented in Fig. |21 The calculated profiles are interpolated 




r(A) 



FIG. 3: Energies of interaction vs. molecule-molecule sep- 
aration for all studied aromatic systems. Calculations are 
performed using BLYP+DCACP functional. 



with cubic splines. The equilibrium separations, r eq and 
the minima of the interaction energy, -E™*, are deter- 
mined from the interpolated curves and are reported in 
TablelU together with the value for two isolated graphcnc 
sheets from Ref. |54j . 

When increasing the disk-size of the systems, the en- 
ergy of interaction per carbon atom increases; corre- 
spondingly, the equilibrium separation decreases. The 
remaining difference of the largest system with respect 
to graphene is most probably due to the static multi- 
pole of the saturating hydrogen atoms, which is known 
to represent up to 7% of the interaction energy in the 
case of benzene [z2- The convergence of the interac- 
tion parameters with the system's size can be exploited 
for extrapolations to even larger structures, such as su- 
pernaphtalene or supertriphenylene [73j |. for which the 
importance of symmetry and hydrogen atoms can be ex- 
pected to decrease, due to the even smaller ratio between 
molecular perimeter and surface. 

To obtain a coarse-grained interaction potential we 
have normalized the energies by the interaction energy 
at the equilibrium distance, -Eg™', and have scaled the 
molecule-molecule separation with the equilibrium dis- 
tance, r eq . After scaling all DFT profiles superimpose on 
a single curve, as illustrated in Fig. 01 This suggests that 
the main contribution to the non-bonded interactions of 
the atoms in the dimer is due to an (additive) pairwise 
carbon-carbon interaction. 

We have fitted the master curve with three frequently 
used potentials, Lennard- Jones (LJ), Morse, and Buck- 
ingham. All the potentials have been constrained to 
have the minimum of the energy, u eq = — 1, at the di- 
mensionless separation, x eq = 1. Under this constraint, 
the Lennard-Jones potential is parameter-free, while the 
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TABLE I: Minimum of the interaction energy, Elf together 
with the equilibrium distance, r eq , for all the systems, graph 
corresponds to the calculated prediction for two graphene 
sheets in vacuum which compares well to experiment. 



System 


N 


-E' e f (kJ/mol) 


r eq (A) 


(A) benzene 


6 


13.8 


3.77 


(B) pyrene 


16 


53.7 


3.67 


(C) triphenylene 


18 


60.6 


3.70 


(D) perylene 


20 


66.3 


3.63 


(E) coronene 


24 


90.1 


3.60 


(F) hexabenzocoronene 


42 


162.7 


3.60 


graph" 


42 


140.7 


3.30 



"from Ref. 154 
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FIG. 4: Interaction energies in units of Elf vs. the separation 
scaled by r eq . The DCACP-BLYP-DFT results (symbols) for 
all systems fall on a single master curve. Solid lines repre- 
sent fits corresponding to different coarse-grained potentials, 
Lennard- Jones, Morse, and Buckingham. 



Morse and Buckingham potentials have each a single fit- 
ting parameter, 



UMorsc(x) 
^Buckingham 



1 

X s ' 
1 - e 
1 

(3-6 



—a(x — l) 



x 6 



(3) 



Here, x — r/r eq . The interval x = [0.9,1.4] has been 
used to determine the fitting parameters a — 5.6 and 
P = 12.3. 

The dimensional potential is obtained by multiplying 
u(x) with the absolute value of E l Jf and using x = r/r eq 
as the argument. The corresponding values are given in 



Table [I] For instance, the Lennard- Jones potential for 
benzene would have a potential energy of interaction of 
13.8 x it lj(?V3. 77), in kJ/mol, and r in A. 

All the three potentials fit very well to the DFT calcu- 
lations. The Morse potential reproduces best the repul- 
sive and the attractive part of the potential. The fact 
that all interaction potentials fall on a single master- 
curve implies that the performance of the fits remains 
constant for all investigated supermolecular systems. 



C. Coarse-grained potentials 

1. United atoms model 

Equations @, can, in principle, be used directly for 
the parametrization of interaction potentials treating the 
whole molecule as one interacting point, such as the Gay- 
Berne potential. Here, however, we will be interested in 
more accurate representations. We start with a united 
atom model, in which all hydrogen atoms are embedded 
into the carbons they saturate. All considered molecules 
are assumed to be rigid, i.e., no stretching, bending, 
or torsional energy is included. Since, already for ben- 
zene, the electrostatic contribution represents up to 7% 
of the total intermolecular energy |x3 , we have neglected 
this contribution for the parametrization of the coarse- 
grained model. 

For atomistic two-body potentials the interaction of 
two molecules is a sum of the corresponding pair inter- 
actions of all atoms. If the effective dispersion-repulsion 
interaction between two atoms (or united atoms) i and j 
of different molecules is represented by a Lennard- Jones 
6-12 potential, 



U ij (r ij )=4e i 



/ \ 12 




&) - 









(4) 



then the molecule-molecule interaction is the sum over 
all pair interactions of the atoms belonging to different 
molecules, 



U(l,2) 



(5) 



In what follows we have assumed that all inner and 
edge carbons have the same parameters, e and a, and 
optimized these parameters to reproduce the desired 
molecule-molecule interaction. For the fit of the DFT 
data to Eq. [51 we have considered only the region close 
to the equilibrium separation r eq . This limitation is due 
to the fact that the employed DCACP-correction to the 
DFT functional was calibrated to reproduce only this 
equilibrium region, and does not explicitly include the 
typical dissociative r _6 -behavior. 

The resulting parameters of the fit are summarized in 
Tab. [HJ The DFT data points, together with the cor- 
responding fitting curves are shown in FiglSJa). Again, 
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FIG. 5: DFT-data and corresponding fitting curves for united 
atoms (a) and benzene-bead (b) representations. The insets 
illustrate the representations for triphenylene (C). 



System 


N 


c 


a 


N B 






(A) benzene 


6 


0.466 


3.556 


1 


13.802 


3.358 


(B) pyrene 


16 


0.408 


3.541 


4 


5.019 


3.428 


(C) triphenylene 


18 


0.407 


3.569 


4 


6.963 


3.417 


(D) perylene 


20 


0.383 


3.512 


4 


4.984 


3.390 


(E) coronene 


24 


0.393 


3.498 


7 


3.536 


3.413 


(F) hbc 


42 


0.353 


3.518 


13 


2.873 


3.459 



TABLE II: Lennard- Jones parameters for united atom (e,a) 
and benzene-bead (es,c"s) parameterizations of all the stud- 
ied systems. The cutoff r cu t = 15 A has been used to evaluate 
the potential in Eq. (|^J. Energies e, es are in kj/mol; a and 
(7b are in A. N is the number of carbons, Nb is the number 
of benzene beads per molecule. 



the Lennard-Jones potential does not reproduce perfectly 
well the attractive tail but, as explained before, the posi- 
tion of the minimum is more important for our purposes. 

We can also compare some of the obtained parame- 
ters to values of existing force fields. For instance, a 
number of united atom force fields are available for ben- 
zene [zS IzM (ZS (Z3, IzM IZ^I ■ Depending on the employed 
parametrization, the literature values are a G [3.25—3.57] 
A and e s [0.4—0.75] kJ/mol, i. e. the values of this study 
fall into the range of parameters predicted from thermo- 
dynamic properties of benzene. 

For the investigated systems, no abundant experimen- 
tal data, particularly concerning dimer interactions, is 
available. However, adsorption energies of benzene or 
coronene on basal planes of graphite were measured |7l| , 
~ 50 kJ/mol and ~ 120 kJ/mol, respectively. While we 
find no quantitative agreement with the results for ben- 
zene, the agreement for coronene is better. It is plausible 
that due to the interaction of a molecule with the bulk 
of graphite, the adsorption energy on graphite is larger 
than the interaction energy between two isolated benzene 



molecules. 



2. Benzene-bead representation 

In coarse-grained simulations, one frequently encoun- 
ters the approximation that fragments building up larger 
systems are rigid, i.e. their internal stretching, bending, 
or torsional energy contributions are neglected. Exploit- 
ing this assumption, a computationally more efficient 
coarse-grained representation can be proposed in which 
each benzene is represented as an interacting point lo- 
cated at its center of mass. This reduces the number of 
degrees of freedom considerably. Assuming, likewise, a 
Lennard-Jones type of interaction, the above presented 
procedure to fit to the DFT-data has been applied to 
obtain the corresponding parameters for benzene-beads. 
The results and fitting curves are displayed in TABLE UTI 
and Fig. [Sf b). respectively. There is no significant dif- 
ference in the interaction profiles upon use of the united 
atom or the benzene-bead representation. However, for 
the latter the size of the beads is relatively small, i.e. only 
as large as the internal diameter of a benzene unit. This 
implies that only mesophases in which molecules experi- 
ence interactions via a face-to-face arrangement can be 
studied. Another limitation is that the rotational profile 
of the interaction energy, e.g. due to azimuthal rotation 
of a moiety, has not been included in the parameteriza- 
tions. Consequently, an accurate prediction of the helix 
structure, often observed in columnar mesophases, can 
not be expected. 



IV. CONCLUSIONS 

The DFT KS-potentials using LDA, BP, PBE, or 
BLYP approximations to the xc-potential fail to correctly 
predict an attractive interaction between polyaromatic 
hydrocarbon molecules. The BLYP-DCACP-extension 
of Ref. [54| has successfully been applied to compute 
interaction energies for hexabenzocoronene, coronene, 
perylene, triphenylene, pyrene, and benzene, without any 
computational overhead or necessity of a priori assign- 
ments of fragments. By scaling the obtained data with 
equilibrium energy and distance values, a single function 
has been found to describe all molecule-molecule interac- 
tions independent of the number of atoms. 

The DFT results have been used to parameterize a 
united atom representation, taking each carbon as an 
interacting site. Additionally, another coarse-grained 
representation has been proposed and parameterized in 
which each benzene unit represents one Lennard-Jones 
bead. 

The obtained interaction potentials will be of use for 
future studies of columnar phases of corresponding com- 
pounds or their derivatives within atomistic molecular 
dynamics simulations or more coarse representations. 
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